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In this paper we explore the practical use of the corner transfer matrix and its higher- dimensional 
generalization, the corner tensor, to develop tensor network algorithms for the classical simulation 
of quantum lattice systems of infinite size. This exploration is done mainly in one and two spatial 
dimensions {Id and 2d). We describe a number of numerical algorithms based on corner matri- 
ces and tensors to approximate different ground state properties of these systems. The proposed 
methods make also use of matrix product operators and projected entangled pair operators, and 
naturally preserve spatial symmetries of the system such as translation invariance. In order to assess 
the validity of our algorithms, we provide preliminary benchmarking calculations for the spin- 1/2 
quantum Ising model in a transverse field in both Id and 2d. Our methods are a plausible alternative 
to other well-established tensor network approaches such as iDMRG and iTEBD in Id, and iPEPS 
and TERG in 2d. The computational complexity of the proposed algorithms is also considered and, 
in 2d, important differences are found depending on the chosen simulation scheme. We also discuss 
further possibilities, such as 3d quantum lattice systems, periodic boundary conditions, and real 
time evolution. This discussion leads us to reinterpret the standard iTEBD and iPEPS algorithms 
in terms of corner transfer matrices and corner tensors. Our paper also offers a perspective on many 
properties of the corner transfer matrix and its higher-dimensional generalizations in the light of 
novel tensor network methods. 
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I. INTRODUCTION 

It is commonly accepted that understanding the prop- 
erties of quantum systems of many particles is one of 
the most important and challenging problems in con- 
densed matter physics. In this sense, much remains yet 
to be understood. For instance, and in spite of many 
efforts, it is still unclear what is the mechanism respon- 
sible for high- Tc superconductivity in cuprates , not to 
mention iron-based and organic superconductors, for 
which this mechanism is also a mistery to a great ex- 
tent. Many other condensed matter phenomena beyond 
Landau's paradigm of phase transitions have also proven 
non-trivial to understand. In this respect, a great in- 
terest in exotic (i.e. beyond Landau's paradigm) phases 
of matter has arised recently. Examples of this are, to 
name a few, topologically ordered phases (where a pat- 
tern of long-range entanglement prevades over the whole 
system)"^, quantum spin liquids (phases of matter that 
do not break any symmetry)^, and deconfined quantum 
criticality (quantum critical points between phases of 
fundament ally- different symmetries) ' . 

The standard approach to understand these systems is 
based on proposing some simplified, physical model that 
is believed to mimic the relevant interactions responsible 
for the observed physics. This is the case of e.g. the 
Hubbard and t — J models for high-Tc cuprate supercon- 
ductors, as well as quantum Heisenberg antiferromagnets 
with frustrating interactions for some magnetic materials 
with a spin liquid ground state . Sometimes we are lucky, 
and these models are exactly solvable. In practice, this 
means that one can compute some (if not all) relevant 
properties analitically. But this is not usually the case 



and, in spite of their apparent simplicity, these models 
turn out to have such an outstandingly complex behavior 
that one needs to rely on alternative approaches. Quan- 
tum simulations, as proposed by Feyman , are certainly 
a possibility. Recent experimental results in this direc- 
tion using e.g. ultracold atoms in optical lattices are in 
fact really promising"^. However, the current technolog- 
ical status does not allow yet to fully understand many 
interesting systems. Thus, one needs to rely on faithful 
methods to implement numerical (classical) simulations. 

From the point of view of numerical simulation al- 
gorithms, there has been increasing interest in recent 
years in the so-called tensor network methods to sim- 
ulate strongly correlated systems . In these methods 
the wave function of the system is described in terms 
of a network of interconnected tensors (a tensor net- 
work). As such, tensor network techniques offer efficient 
descriptions of quantum many-body states of the system 
that are based on the amount of entanglement in the 
wave function. The amount and structure of entangle- 
ment is a consequence of the chosen network pattern and 
the number of parameters in the tensors. The most fa- 
mous example of a tensor network method is probably the 
Density Matrix Renormalization Group (DMRG) , 
introduced by White in 1992, and which has been the 
technique of reference for the last 20 years to compute 
ground state properties of Id quantum lattice systems. 
Recently, though, many important insights coming from 
quantum information science have motivated the appear- 
ance of a host of new tensor network methods. Now- 
days it is easy to get lost in the zoo of names for all 
these methods, e.g. Time-Evolving Block Decimation 
(TEBD)'^ Folding Algorithms'^, Projected Entagled 
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Pair States (PEPS)^^, Tensor-Entanglement Renormal- 
ization Group (TERG)^^, Tensor Product Variational 
Approach , Weighted Graph States^^, Entanglement 
Renormalization (ER) ' , String-Bond States ' , 
Entangled-Plaquette States , Monte Carlo Matrix Prod- 
uct States^'', Tree Tensor Networks^^"'^^ and Continuous 
Matrix Product States , just to name some of the most 
recent ones. Each particular method has its own ad- 
vantages and disadvantages, as well as optimal range of 
applicability. 

Tensor network methods are also an interesting ap- 
proach to many-body systems since they offer a lot 
of flexibility. For instance, with tensor networks one 
can study a variety of systems in different dimensions, 
of finite or infinite size with different boundary 

conditions"^^'^^ , symmetries^^"'^^, as well as systems of 
bosons^ fermions and frustrated spins^^"^^. Dif- 
ferent types of phase transitions as well as the robust- 
ness of topological order to local perturbations^^"^"^ have 
also been studied in this context. The fact that it is possi- 
ble to develop algorithms for infinite-size systems is quite 
relevant, since this means that it is possible to approxi- 
mate the properties of a given system directly in the ther- 
modynamic limit and without having to rely on finite-size 
scaling extrapolations. This is achieved by cleverly ex- 
ploiting the translational invariance of the system. Ex- 
amples of methods using this approach are iDMRG'^'^'^"^ 
and iTEBD in Id (the 'i' means infinite), as well as 
iPEPS^^'^o and TERG^'^ in 2d. 

The possible variations of all these methods are, in 
practice, unending. In this respect, we feel that a use- 
ful tool still relatively unexplored in this context (albeit 
with exceptions) is the so-called Corner Transfer Matrix 
(CTM)^^'''^. This was originally introduced by Baxter 
in 1968 in the context of classical statistical mechanics 
and exactly solvable models ' . The CTM has very 
nice properties, specially regarding its spectrum of eigen- 
values, and has been a standard tool to find the exact 
solution of many classical 2d models such as the hard- 
hexagon and related models ' . However, its practical 
use goes well beyond analytical solutions, and numerical 
algorithms can be developed to approximate the proper- 
ties of 2d classical lattice models based on clever trunca- 
tions in the eigenvalue spectrum of CTMs. Baxter him- 
self was the first to explore this possibility by means of 
a variational method '^'^^ that was an extension of the 
so-called Krammers-Wannier approximation^ . In fact, 
Baxter's method can be understood as a precursor of 
DMRG but in the context of classical lattice models (this 
statement will be made more precise later on). The for- 
mal combination of CTMs and DMRG was later on put 
forward by Nishino and Okunishi in the so-called Corner 
Transfer Matrix Renormalzation Group (CTMRG) ' . 
Since then, numerical algorithms using CTMs have been 
widely used, mostly focusing on classical lattice mod- 
els. From the point of view of quantum lattice systems, 
Ref. ' ' already discussed the possibility of simulating Id 
quantum systems by using CTMs and a Suzuki- Trotter 



decomposition of the evolution operator^^'^^. Neverthe- 
less, we believe that CTMs have not yet been fully ex- 
ploited as a tool in the context of the novel tensor network 
methods that are being developed for quantum lattice 
systems. 

Our aim in this paper is to cover in part this gap by 
exploring the applicability of the CTM to develop algo- 
rithms for the simulation of quantum lattice systems. We 
do this mainly in Id and 2d, which naturally leads us to 
consider the generalization of the CTM to higher dimen- 
sional systems. Following the convention from previous 
works ' ' , we call this generalization corner tensor. 
More specifically, in this paper we describe a number of 
numerical algorithms based on CTMs and corner ten- 
sors to approximate ground state properties of quantum 
lattice systems of infinite size. The methods that we 
present rely also on the use of Matrix Product Operators 
in Id and Projected Entangled Pair Operators in 2d, and 
naturally preserve the spatial symmetries of the system, 
including invariance under translations. This, in turn, 
is a significant difference with respect to some previous 
approaches for infinite systems that slightly break trans- 
lational invariance Also, we will see that these 
methods produce, in a natural way, individual tensors 
for the 'bra' and 'kef parts of local expectation values, 
in a way similar to the so-called single-layer picture^^. 
This does not seem to have remarkable consequences in 
Id but, as we shall discuss, it has some interesting im- 
plications for the 2d algorithms when compared to other 
methods such as e.g. iPEPS ' . 

In order to prove the validity of our algorithms we pro- 
vide preliminary benchmarking calculations for the spin- 
1/2 quantum Ising model in Id and 2d. Our methods are 
roughly comparable in accuracy to other well-established 
approaches such as the iDMRG ' and iTEBD "^'^^^^ 
methods in Id and the iPEPS method and TERG in 
2^19,37,70^ thus offering a possible alternative to them. 
Moreover, we will also discuss that it is possible in prin- 
ciple to use some of the 2d results for the development 
of a 3d algorithm in combination with some tensor up- 
dates, in the same spirit as the algorithm in Ref.^^ for 2d 
systems. The computational complexity of all the pro- 
posed algorithms is also analyzed and, as we shall see, 
important differences are found depending on the dimen- 
sionality and chosen simulation scheme. 

For completeness, we also discuss briefly in Appendix 
B the case of periodic boundary conditions, as well as real 
time evolution. This discussion will lead us to a beautiful 
interpretation of the standard iTEBD and iPEPS algo- 
rithms in terms of CTMs and corner tensors. We also 
believe that the algorithms described in this paper will 
be useful for the pracical implementation and develop- 
ment of further tensor network methods in the future. 

This work is organised as follows: In Sec. II we present 
several preliminary concepts and generalities. These in- 
clude notions on tensor networks (Sec. II. A), corner trans- 
fer matrices (Sec.II.B,C), Matrix Product States and 
Projected Entangled Pair States (Sec. II. D), and time evo- 
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FIG. 1: (color online) (a) Tensor network diagrams repre- 
senting a scalar, a vector, a matrix and a three-index tensor, 
(b) The grouping of two indices a and /3 of a tensor pro- 
duces another index represented as a thick line, (c) 
Diagram representing a matrix multiplication, or contraction 
of an index, (c) Diagram for tensor network A/", with four 
open (non-contracted) indices a, /3, 7 and b. The result of the 
contraction is a four-index tensor T. 



lution with Matrix Product Operators and Projected En- 
tangled Pair Operators (Sec.II.E). In Sec. Ill we present 
two simple algorithms based on CTMs and corner ten- 
sors to approximate ground state properties of infinite- 
size quantum lattice systems. More elaborated versions 
of these algorithms are presented in Appendix A. After 
discussing the general approach in Sec. III. A, we consider 
\d systems in Sec. III. B and 2d systems in Sec. III. C. A 
summary of the proposed methods and their complexi- 
ties is done in Sec. III. E. In Sec. IV we present preliminary 
benchmarking numerical calculations for \d and 2d sys- 
tems. Sec. V contains our conclusions and final remarks. 
Moreover, in Appendix B we discuss further possibilities, 
such as Zd quantum lattice systems, periodic boundary 
conditions, and real time evolution. 



II. PRELIMINARY CONCEPTS 

The goal of this section is to introduce some prelim- 
inary concepts that we will need in the presentation of 
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FIG. 2: (color online) Different orderings in the contraction 
lead to different number of operations. Assuming that all the 
indices of the tensors can take up to x different values, then 
the number of operations is scheme (a) is O(x^), whereas it is 
O(x^) in scheme (b). Thus, scheme (b) is more efficient than 
scheme (a). 



the algorithms of the forthcoming sections, as well as to 
provide some background on relevant topics. We start 
by quickly reminding some basic concepts on tensor net- 
works and tensor network diagrams. Then, we review 
the basics of the CTM and its properties, as well as sev- 
eral aspects regarding possible gneralisations. After this, 
we quickly remind the basic definition of Matrix Prod- 
uct States (MPS) and Projected Entangled Pair States 
(PEPS), and then review the implementation from Ref.^-^ 
of time evolution operators (both in real and imaginary 
time) using Matrix Product Operators (MPO) and Pro- 
jected Entangled Pair Operators (PEPO). 



A. Tensor networks and diagrams 

For our purposes, a tensor is a multidimensional ar- 
ray of complex numbers. A Tensor Network (TN) is a 
network of tensors whose indices are connected according 
to some pattern. This connection of indices is done by 
summing over all the possible values of common indices 
between tensors. Summing over an index is also called 
contracting the index. Summing over all the possible in- 
dices of a given TN is called contracting the TN. 

Instead of using equations, tensors and TNs are more 
easily handled by using a diagrammatic notation in terms 
of tensor network diagrams^ see Fig.(l). In these dia- 
grams tensors are represented by shapes, and indices in 
the tensors are represented by lines emerging from the 
shapes, see Fig.(l.a,b). A TN is thus represented by a 
set of shapes interconnected by lines. The lines connect- 
ing tensors between each other correspond to contracted 
indices, whereas lines that do not go from one tensor 
to another correspond to open indices in the TN, see 
Fig.(l.c,d). As expected, the contraction of a TN with 
some open indices gives as a result another tensor, and 
in the case of not having any open indices the result is 
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FIG. 3: (color online) Partition function Z(I3) at inverse tem- 
perature /3 of some classical lattice model on the square lat- 
tice, defined by a symetric weight matrix T between nearest 
neighbours (such as e.g. the classical Ising and Potts mod- 
els). The partition function can be written as a TN with 
either one tensor T per link, or one tensor per site resulting 
from the contraction of four VT tensors (see e.g. Ref. ' ) 
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FIG. 4: (color online) (a) Matrix Product State (MPS), (b) 
Matrix Product Operator (MPO). (c) Projected Entangled 
Pair State (PEPS), (d) Projected Entangled Pair Operator 
(PEPO). The physical dimension is q in all cases, whereas 
the bond dimension is x for MPS and MPO, and D for PEPS 
and PEPO. 



a scalar. Notice, though, that the total number of op- 
erations that must be done in order to obtain the final 
result of the contraction depends strongly on the order 
in which tensors in the TN are contracted, see Fig. (2). 
To minimize the computational cost of a TN contrac- 
tion, one must thus optimize over the different possible 
contraction orderings. 

There are famous examples of TN in the context of 
many-body physics. For instance, the partition function 
of a (i- dimensional classical lattice model with nearest- 
neighbour interactions is a TN in d dimensions (see 
Fig. (3)). Also, for quantum lattice systems, the classes of 
MPS and PEPS are suitable to describe quantum states 
of Id and 2d systems respectively (see Fig. (4. a, c)). Other 
examples make use of an extra 'holographic' dimen- 
sion accounting for some renormalization group scale^^, 
such as Tree Tensor Networks (TTN) (Fig.(5.a)) and 
the so-called Multi Scale Entanglement Renormalization 
Ansatz (MERA) (Fig.(5.b)), which is at the basis of En- 
tanglement Renormalization'^'^''^'^. TNs can also be used 




FIG. 5: (color onhne) (a) Tree Tensor Network (TTN). (b) 
Id binary Multi Scale Entanglement Renormalization Ansatz 
(MERA), with unitaries u and isommetries w. For details, 
see Ref.22'23,28-3 1 ^ 



to describe operators, such as MPOs in Id and PEPOs 
in 2d (Fig.(4.b,d)). 

B. Corner Transfer Matrix: fundamentals 

The concept of Corner Transfer Matrix (CTM) was 
introduced by Baxter in the context of exactly solvable 
models in 2d classical statistical mechanics ' . These 
can be defined for any planar TN. However, for the sake 
of simplicity, here we shall assume the case of a 2d TN 
on a square lattice as in Fig. (6). This TN could repre- 
sent, for instance, the partition function of some 2d clas- 
sical lattice model as in Fig. (3) or, as we will consider 
in Sec. Ill, the imaginary time evolution of a quantum 
Id system. In order to define the CTM one makes the 
following observation: the contraction of the TN can be 
obtained by multplying four matrices Ci, C2, C3 and C4, 
one for each corner (see Fig. (5. a)). Thus, 

Z = tr(CiC2C3C4) , (1) 

where Z is the scalar resulting from the contraction. 

Matrices Ci , C2 , C3 and C4 are the Corner Transfer 
Matrices of the system. They correspond to the contrac- 
tion of all the tensors in each one of the four corners of 
the 2d lattice of tensors. The nomenclature 'transfer ma- 
trix' is appropriate, since the CTM 'transfers' a vector 
in the angular direction around the center of the lattice 
by an angle of 7r/2 in our case. To further simplify our 
explanation, let us assume that the four CTMs are equal, 
that is C = Ci = C2 = C3 = C4. 

Sometimes it is convenient to define diagonal CTMs 
Cd = PCP~^. Depending on the symmetries of the sys- 
tem (and thus of C), matrix P may be arbitrary, unitary 
or orthogonal. Let us assume that there are x different 
eigenvalues with a = 1,2,.. .,x. Then, the contrac- 
tion of the full TN adopts the very simple expression 

Z = tr{Cj) = j2'^t- (2) 




FIG. 6: (color online) (a) The contraction of a 2d square 
lattice of tensors results in a scalar Z. This contraction can 
be understood as the trace of the product of four CTMs, one 
for each corner. Notice that the exact CTMs will, in principle, 
have indices that can take up to exponentially many different 
values (in the size of the system). Thus, their indices are 
represented by a thick black line. The aim of CTM numerical 
methods is precisely to reduce the size of this index in some 
optimal or quasi-optimal way. (b) A possible reduced density 
matrix p of a system with the same CTM at every corner. 
Open indices of p are shown in red. 



The fact that we choose to name the number of different 
eigenvalues as x is made on purpose. There is in fact a 
direct relation between this parameter and the so-called 
bond dimension of an MPS. The relation between these 
two seemingly different parameters will be made clearer 
later on in Sec.IILB, when considering Id quantum lat- 
tice systems. 

CTMs are of paramount importance in the context 
of classical statistical mechanics. They have been used 
to solve the hard hexagon model, as well as many 
others^^'^"^. But they have also been useful in the context 
of quantum information, since it is known from long ago 
that their eigenvalue spectrum can be related to the en- 
tanglement spectrum of Id quantum lattice systems, see 
e.g. Ref. ' . As we will see, this is a key property in 
the algorithms that will be further described in SecIITB. 
Also, from a numerical perspective, a variational method 
to approximate the partition function per site of a 2d 
classical lattice model was developed by Baxter by trun- 
cating in the eigenvalue spectrum of the CTM ' . This 
idea was later on refined by Nishino and Okunishi, who 
developed the so-called Corner Transfer Matrix Renor- 
malization Group (CTMRG) ' . As such, CTMRG 
is an algorithm to approximate properties of 2d classi- 
cal lattice models with isotropic interactions (and thus 
a high degree of symmetry), and runs by truncating in 
the largest eigenvalues in magnitude of the spectrum of 
matrix C^, which is a reduced density matrix of the sys- 
tem (see Fig.(6.b)). This spectrum, in turn, is given by 
the numbers in Eq.(2), and is thus in one-to-one cor- 
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FIG. 7: (color online) Same as in Fig. (6), but for a cubic 
lattice in 3d. (a) A 3d cubic llatice of tensors results in a scalar 
Z. This contraction can be understood as the contraction of 
eight corner tensors. In the exact case, the corner tensors 
have 'fattened' indices, exactly as in the 2d case. We have 
chosen not to draw explicitely the 3d cubic lattice in order to 
keep the figure simple, (b) A possible reduced density matrix 
p of a system with the same corner tensor at every corner. 
Open indices of p are shown in red. 



respondence with the spectrum of the CTM C. 

As we shall explain later on, from the point of view 
of quantum states of Id quantum lattice systems the 
numbers are, in fact, the spectrum of eigevalues of 
the reduced density matrix of half an infinite chain. Or 
what is the same, the spectrum of Schmidt coefficients 
of half an infinite chain is given by = i^a- Thus, in a 
way, we can think of Baxter's variational method as sort 
of a precursor of the truncation scheme used in DMRG 
and TEBD, but in the context of 2d classical statistical 
mechanics^ Also, it is well known that these CTM 
methods work assymptotically in the limit of a system of 
infinite size. In fact, the closer the method is to conver- 
gence, the more faithful are the truncations in the eigen- 
values of the corresponding reduced density matrix. The 
convergence and speed of CTMRG for classical 2d mod- 
els is also remarkably better than other approaches such 
as the so-called Transfer Matrix Renormalization Group 
(TMRG)^^'^°^'i°^ This is understandable, because in 
many situations away from criticality the largest eigen- 
value of the CTM is non-degenerate, and there is a rather 
big gap to the next lower eigenvalues (as compared to the 
gap of the row-to-row transfer matrix) This ensures 
a fast numerical convergence of CTM methods. 
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C. Corner Transfer Matrix: generalizations 

Since the foundational works by Baxter, there have 
been several attempts to generalize the CTM in a num- 
ber of different ways. Here we mention some basic pos- 
sibilities. 

The first obvious generalization is for 2d lattices dif- 
ferent from the square one. Of course, the CTM can be 
defined for any planar lattice as long as corners can be 
defined. This includes the usual lattices, but also more 
exotic constructions such as lattices with hyperbolic ge- 
ometry (or negative curvature), e.g. lattice discretiza- 
tions of Anti de Sitter (AdS) manifolds . 

An attempt to improve the efficiency of the truncations 
involved in CTMRG was also proposed in Ref. , where 
a directional version of the method was put forward in 
the context of 2d iPEPS calculations, and where tensors 
were no longer real and broke rotational symmetry. In 
Sec. III. B of this paper we will also propose another plau- 
sible way of doing this generalization to non-symmetric 
tensors in the context of the algorithms that we explain 
here. 

Another natural generalization is the case of periodic 
boundary conditions. Even if this sounds counterintu- 
itive (since a periodic system does not have any corners) , 
it turns out that the ideas from CTMRG can be general- 
ized to deal with this type of systems as well . We will 
briefly comment on how one can perform this generaliza- 
tion in Sec. III. D. 

Moreover, the CTMRG technique can also be general- 
ized to deal with stochastic models, see Ref. for details. 
Also, a quantum counterpart of the CTM is the so-called 
Corner Hamiltonian, which essentially is the logarithm of 
the CTM. Numerical techniques for Corner Hamiltonians 
have also been developed, see e.g. Ref. . 

Finally, the CTM can also be generalized to higher- 
dimensional systems. The natural corner object now is 
not a matrix, but a tensor with three indices (for a cubic 
lattice) which we call the corner tensor , see Fig. (7). 
In terms of corner tensors Ci, C2, . . . , Cg, the contraction 
of a 3d lattice of tensors reads 



^ — / (C'l, C2, C3, C4, C5, Ce, C7, Cg) , (3) 

where / is a function that performs the corresponding 
TN contraction (see Fig. (7. a)). This higher-dimensional 
generalization is quite obvious but, yet, it has some non- 
trivial consequences since one is dealing now with tensors 
instead of matrices. In particular, the corner tensors can 
no longer be diagonalized (or more precisely, there is no 
unique eigenvalue decomposition since it depends on how 
one chooses to split the indices of the tensor). Thus, no 
expression like the one in Eq.(2) for the 2d case can be 
obtained in general for 3d. Nevertheless, corner tensors 
still have interesting spectral properties with respect to 
eigenvalue/singular value decomposition of bipartitions 
of their indices. The behavior of these singular values will 



be the key to define simplified numerical approaches for 
higher-dimensional systems, as we shall do in Sec. III. C. 

D. Matrix Product States and Projected 
Entangled Pair States 

Let us now revise briefly the concepts of Matrix Prod- 
uct States (MPS) and Projected Entangled Pair States 
(PEPS). There is a vast amount of literature on these 
two families of states, and we refer the interested reader 
to it for further details (see e.g. Ref. and references 
therein). 

Consider a quantum many-body system of N particles. 
The quantum state of the system is |^) G where H = 
^^i'^^^^ is the total Hilbert space of the system and 
H^^] is the individual Hilbert space of particle r. Let us 
assume that each particle is modelled by a g-level system, 
so that dim(H[^^) = q. Given a local basis |v) for each 
site r, with v = 1,2,...,^, the quantum state of the 
system reads 

1^) = ^iii2-iN\hi2---iN)' (4) 

iii2---iN 

The coefficient Ci^i^...i^ can be understood as a tensor of 
N indices with 0{q^) complex coefficients. Thus, this 
is clearly an inneficient representation of the quantum 
state because the required number of parameters scales 
exponentially with the size of the system. In order to find 
an efficient description, one can consider a decomposition 
of the above tensor into a MPS for Id, or a PEPS in 2d. 
These decompositions are shown in the tensor network 
diagrams of Fig. (4. a. c), where open boundary conditions 
are assumed. 

Both MPS and PEPS offer an efficient description of 
the quantum state |^) for Id and 2d respectively. More- 
over, MPS and PEPS are known to have many interesting 
(and important) properties. For instance, both of them 
satisfy the so-called area law scaling of the entanglement 
entropy^^"^^, which is a key property of low-energy states 
of most quantum many-body systems (albeit with no- 
table exceptions, see Refs. ' ). Moreover, it is known 
that ground states of Id gapped local Hamiltonians can 
be efficiently approximated with exponential accuracy by 
an MPS , and the same holds for thermal states in 2d 
with PEPS . From the numerical point of view, MPS is 
the relevant class of states at the heart of efficient meth- 
ods for Id systems such as DMRG and TEBD. PEPS 
is also at the basis of simulation methods for 2d systems 
such as the finite and infinite PEPS methods ' , as well 
as TERG ^ 

We wish to remark here a couple of properties of MPS 
and PEPS. First, for a system of size the number 
of parameters in both families of states scales linearly 
with and polynomially in the bond dimension of the 
tensors (that is, the number of different values for the 
connecting indices in the MPS or PEPS tensor network, 
see Fig. (4)). We call this bond dimension x in the case of 
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MPS, and D in the case of PEPS. Importantly, this bond 
dimension can be regarded as a measure of the number of 
parameters in the TN, but also as a measure of the maxi- 
mum amount of entanglement that can be handled by the 
wavefunction (see Ref. ). Second, for systems invariant 
under translations, it is possible to take the thermody- 
namic limit N ^ oo and consider a MPS or PEPS for 
a system of infinite size. This is done by repeating the 
same unit cell of tensors across the whole lattice'^^"^^. 
This trick is at the basis of methods to study infinite- 
size systems such as iDMRG, iTEBD and iPEPS. The 
methods that we shall propose in this paper will be for 
infinite-size systems, and thus rely on this idea as well. 

E. Time evolution with Matrix Product Operators 
and Projected Entangled Pair Operators 

Let us now consider the problem of time evolution. 
We assume that this evolution is generated by a Hamil- 
tonian H that is the sum of local interaction terms on lat- 
tice. For simplicity of the explanation, let us also assume 
a time-independent Hamiltonian with nearest-neighbour 
interactions, namely 

(r,r') 

although time-dependent cases and more generic types 
of intereactions could also be considered (including long- 
range ones ' ). The real time evolution of a given state 
1^(0)) reads 

\^{t)) = e-^^'\m)- (6) 

It is equally possible to consider the evolution in imagi- 
nary time r in order to get the ground state |^^s) of 
namely 



Ref.^^, but which are important for our purposes. For 
concreteness of the explanation, let us imagine that H 
corresponds to the spin- 1/2 ferromagnetic quantum Ising 
model in a homogeneous transverse field, 

(r,r') r 

where and are the usual Pauli matrices, and the 
sum in the first term is over nearest-neighbours. Let us 
now assume, for simplicity, the case of a Id system of 
N particles with periodic boundary conditions. With- 
out loss of generality, we also consider the evolution in 
imaginary time under this Hamiltonian for a total time 
T. The first thing to do, and as also done also in other 
approaches^^"^''''''^, is to approximate the whole evolution 
by breaking the total evolution time T into smaller steps 
of size (^r <C 1. Then, we have that the evolution oper- 
ator U = exp(— i^T) can be written as /7 = [U{Sr)]'^ = 
[exp(— i^(5r)]^, with m = T/Sr. Next, we would like to 
describe the term U{Sr) = exp(— i^^r) by an MPO. We 
wish to make it in such a way that the resultant MPO is 
invariant under translations and symmetric under space 
inversion, since these are symmetries of the original 
Hamiltonian, and also favors stabiltity in numerical ma- 
nipulations. As explained in Ref. , this can be achieved 
in the following way: first, split H = + Hx, so that 
we can perform a Suzuki- Trotter decomsposition^^'^^ 
exp(-i^^r) ^eM-HJr)exp{-H,dr)^0{6T^y^\ In 
this splitting, contains all the terms with oper- 
ators and Hx the terms with ax operators. The term 
ex.-p{—HxSr) is just a tensor product of one-body opera- 
tors, 

and its MPO representation is trivial (see Fig. (8. a)). 
However, the term ex.p{—HzSr) requires a little bit more 
consideration. In order to evaluate the MPO for this 
term, we first remember that 



^^gs) = lim 



e-^"|^(Q)) 
||e-^1^(0))|| ' 



(7) 



where we have assumed that the initial state |^(0)) has 
a non-zero overlap with the ground state |^^s), and the 
appropriate normalization of the state has been included. 

As is well known, these two different types of evolutions 
can be approximated using a TN approach for Id and 2d 
systems, both for finite and infinite systems. What is 
important for us, though, is that the evolution operators 
e~*^^ and e~^^ can usually be approximated by a se- 
quence of well-behaved MPOs in Id or PEPOs in 2d^ 
as explained in Ref. . This MPO and PEPO approxi- 
mation will turn out to be a key step in the algorithms 
that will be proposed in Sec. HI. 

Let us now remind the basic steps in obtaining MPO 
and PEPO descriptions for the evolutions generated by 
H. Here we simply review some of the results from 



= COSh {uj)lW'^'^ + siuh {uj)aPa^^'^ (10) 

for any u. Using this property, and a little bit of algebra, 
it is easy to arrive at an expression for the MPO. If the 
operator is expressed as 

i' s,j' s 

(11) 

where \ir) is the basis of eigenstates of at site r (and 
same for the 7's), then the coefficients c^/\^"'Y are given 
by the MPO 

where the non-zero components (Mj)^^ of tensor M are 
(Mj)n = cosh (St)!), {M})22 = sinh {5t)I), 
(Mj)i2 = (Mj)2i = Vsinh((5T)cosh((5T)(<j,)i , (13) 
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FIG. 8: (color online) (a) Diagrammatic representation of the 
tensor M in Eqs. (12,13). Bond indices are highlighted in pink, 
(b) Construction of the Id MPO. The contraction of tensors L 
and M results in the MPO tensor R. (c) Construction of the 
2d PEPO. The contraction of one tensor L and two tensors M 
(one along each direction of the lattice) results in the PEPO 
tensor R. 




(b) 




St 



IV^(O))' 
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FIG. 9: (color online) (a) Time evolution of an MPS for |^(0)) 
as driven by an MPO for U{6t). (b) Time evolution of a PEPS 
for 1^(0)) as driven by a PEPO for U{St). 



as explained in Ref. . This MPO is represented in 
Fig.(8.b). Then, in order to obtain a complete MPO 
combining the and terms, we just need to com- 
bine the tensors for the two pieces as shown in Fig.(8.b). 
This MPO description is particularly convenient for nu- 
merics, since it is symmetric with respect to space inver- 
sion, the tensors are real numbers, and it is also invariant 
under translations. Notice also that one can easily take 
the thermodynamic limit N ^ 00. 

It is a good idea to make the tensors in the MPO de- 
scription real and as symmetric as possible. This can usu- 
ally be achieved by a variety of tricks depending on the 
system considered. For instance, as explained in Ref. , 
for the antiferromagnetic Heisenberg models on bipartite 
lattices it may be convenient to perform a sublattice ro- 
tation prior to finding the MPO representation in order 
to make the MPO tensors real. It is a good idea to ap- 
ply tricks like this also to other models, whenever this is 
possible. 

In 2d it is equally simple to find a PEPO descrip- 
tion for the evolution operator. The construction is 
based on that of the MPO for Id. Let us imag- 
ine that this time we have the same Hamiltonian as 
in Eq.(8) but in a square lattice. Then, we can de- 
compose the evolution term generated by Hz into a 



sum of two contributions: one for interactions along 
the rows, and one for interactions along the columns, 
H = i^^^^ + H'^^'K Then, we have that exp(-i^^^r) = 
exp(— iJJ^'^^r) exp(— i^^^^(5r), where 77^^^/^^^ contains 
the interactions alog the rows/columns. It is clear that, 
for each individual row or column, one can find MPO 
descriptions as the ones described before for the Id case. 
Therefore, it is just a matter of combining together all 
these MPOs in a appropriate way in order to obtain a 
PEPO for the desired evolution. This is represented in 
the diagram of Fig. (8. c). This PEPO has again nice prop- 
erties: it is real, and also is symmetric with respect to 
space inversions in the two lattice directions. 

The construction of MPOs and PEPOs that we have 
reviewed here can be further generalized to a variety 
of other models and interactions. Notice that it is of 
course possible to use other alternative approaches to 
build MPOs and PEPOs for evolution operators. How- 
ever, these constructions may not always guarantee the 
symmetry conditions of the obtained TN such as space in- 
version and translational invariance. As we shall see, it is 
important for our purposes that the obtained MPOs and 
PEPOs have these nice symmetries, since this improves 
the efficiency of the algorithms that we shall propose in 
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Sec.IIL 



III. ALGORITHMS 

In this section we propose a number of numerical sim- 
ulation algorithms that allow to compute ground state 
properties of quantum lattice systems of infinite size. Our 
main focus is on Id and 2(i, but we also discuss briefiy 
the possibility of 3d systems, periodic boundary condi- 
tions, and real time evolution. Our algorithms are valid 
to deal with systems with a different amount of spatial 
symmetry. As we shall see, more spatial symmetry means 
simpler and more efficient simulation methods. In Id this 
may not be too relevant, since the number of operations 
in the proposed methods differ only in subleading and 
multiplicative constant terms. But in 2d this is crucial 
since the difference in the leading number of operations 
turns out to be huge, in fact several orders of magnitude, 
depending on the amount of symmetry. The algorithms 
that we present here allow to compute ground state prop- 
erties. As we shall discuss more quantitatively in Sec. IV, 
these methods are a possible alternative to other meth- 
ods for systems of infinite size such as iDMRG, iTEBD, 
iPEPS and TERG^^'^^-^^'^°. Importantly, ah the algo- 
rithms of this section preserve at every step the invari- 
ance under translations. 

First we discuss the general approach, and then we 
describe the details of each algorithm. The cases of 3d 
and periodic boundary conditions are briefiy discussed 
in Sec. III. D. For an overall view of the complexity of the 
different the methods, one can jump directly to Sec. III. E. 



A. General approach 

The general idea is quite simple, and can be under- 
stood from the diagrams in Fig. (9) and Fig. (10). The 
goal is to compute the evolution of some quantum state 
as driven by some Hamiltonian with local interactions 
for a total time T. As we have explained in the pre- 
vious section the time evolution operator, both in real 
and imaginary time, can be decomposed as the action of 
MPOs or PEPOs for a given time interval St on a given 
quantum state described by an MPS or PEPS. Thus, the 
whole evolution of the system can be represented by some 
initial MPS or PEPS, to which one applies the MPO or 
PEPO driving the evolution for as many steps m = T/5r 
as needed, see Fig. (9). In the case of a ground state cal- 
culation the number of steps m is infinite or, in practice, 
very large until convergence of some relevant quantity is 
achieved. 

The main goal of our algorithms is the efficient approx- 
imation of expectation values of local obsevables. For 
instance, let us consider a one-body operator oi. Let us 
also consider the evolved state |^(T)) at some given time 
T, as represented in Fig. (10) for an MPS (the case of a 
PEPS is a straightforward generalization). The expec- 




FIG. 10: (color online) Diagram for the numerator in Eq.(14). 
The expression is hermitian with respect to the vertical in- 
dices, hence the TN is symmetric with respect to a specular 
reflection across the red dashed line plus complex conjuga- 
tion. This means that the individual tensors are symmetric 
with respect to transposition of the vertical indices plus com- 
plex conjugation. 



tation value of operator oi in the evolved state is given 
by 

^ (^(T)loil^(T)) ^ (^(0)|[/tox?7|^(0)) 
^ (^(T)I^(T)) (^(0)|/7t/7|^(0)) ' ^ ^ 

where U is the corresponding evolution operator that can 
be decomposed as = [U{6r)]^^ with U{Sr) being ap- 
proximated by an MPO or PEPO. Thus, e.g. the nu- 
merator in Eq.(14) can be represented diagramatically 
as in Fig. (10). The actual expectation value is the ra- 
tio of two contracted TNs like the one in the figure, one 
with the obserble Oi for the numerator, and one with- 
out the observable for the denominator. For a quantum 
lattice system in d spatial dimensions, these are {d -\- 1)- 
dimensional TNs, where the extra vertical dimension is 
time. Therefore, for an MPS one has to deal with the 
contraction of a (1 + l)(i TN, whereas in 2d one has the 
contraction of a (2 + l)d TN. 

All the existing TN methods that compute real and 
imaginary-time evolutions approximate, in one way or 
another, a contraction like the one described above or 
similar. For instance, iTEBD and iPEPS methods com- 
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(a) 




hermitian 



FIG. 11: (color online) (a) Symmetries of the MPO tensor 
R for the simplfied one-directional Id method. In the gen- 
eral case (full one- and two-directional Id methods), only the 
first equality holds, (b) TN structure for the simplified one- 
directional Id method, (c) General TN structure for the full 
one- and two-directional Id methods. 

pute an MPS or PEPS approximation to the dominant 
eigenvector of some transfer matrix operator defined as 
the MPO or PEPO driving the evolution^^-^^'^^. But 
there are more ways of contracting these TNs. For in- 
stance, one could think of 'folding' the TN across the 
red dashed line in Fig. (10) accompanied by a transversal 
contraction, as has been done in Id in Ref. . 

Our approach in this paper is to use CTMs and corner 
tensors to approximately contract TNs like the one in 
Fig. (10) for different types of systems. This is actually a 
direct generalization of the ideas of CTMRG in Ref.^^'^^. 
In our case, though, it is worth stressing a peculiarity 
inherent to the case of quantum lattice systems: any ex- 
pectation value is by construction a hermitian expresison. 
This means that the coresponding TN is symmetric with 
respect to transposition of the vertical indices and com- 
plex conjugation^ see Fig. (10). Even if quite obvious, this 
property needs to be explicitely taken into account in all 
our algorithms. 

The methods that we propose are different ways of ap- 
proximating these expectation values by using CTMs and 
corner tensors. This is done in two steps. First, one finds 
a set of renormalized tensors (e.g. renormalized CTMs) 
accounting effectivly for the most important correlations 
in the TN. Second, expectation values are approximated 
by using these renormalized tensors. In what follows we 
show how to do this in Id and 2d. 



B. Id quantum lattice systems 

In this paper we consider three algorithms using CTMs 
to approximate expectation values as in Fig. (10). The 
first algorithm is called 'simplified one-directional Id 
method', and is very efficient. This approach is de- 



signed for systems such that the MPO is also hermitian 
with respect to the horizontal indices and is explained 
here. The other two algorithms are considered in Ap- 
pendix A.a,b. The second algortithm is the 'full one- 
directional Id method', valid for an MPO without the 
previous symmetry requirement. The third algorithm is 
the 'full two-directional Id method', which is equivalent 
to the CTMRG algorithm generalized to deal with the 
MPO considered here without extra symmetry require- 
ments. As we shall see, the notion of 'one-directional' 
or 'two-directional' for each approach refers to the num- 
ber of directions in which the lattice is simultaneously 
expanded at every step. 

The relevant parameter in these methods will be 
which is the rank of the renormalized CTMs at every 
step. From Eqs.(l,2) and Figs.(12.d,13.d,25.c,26.c,28), 
one can see that this is also the rank of the the renor- 
malized reduced density matrices of the system. Thus, 
some of the truncations performed by the algorithms in 
this section are really truncations in the entanglement 
spectrum (or, equivalently, the spectrum of Schmidt co- 
efficients) of half an infinite chain 

The leading scaling of the running time (or complexity) 
of the three methods is the same, namely 0{x^). This 
is exactly the same complexity as iDMRG and iTEBD 
methods. However, the multiplicative corrections to this 
overall complexity are very different for each case, which 
implies that the total number of operations is also quite 
different. The simplfied one-directional Id method is the 
less time-consuming, whereas the full two-directional Id 
method is the one that takes more time. 

The procedure in these algorithms is similar to the 
directional CTM algorithm from Ref. . Namely, we (i) 
insert rows and columns in the TN in order to expand its 
structure, (ii) absorb the inserted tensors by contractions 
towards the horizontal {x) and vertical {y) directions, or 
towards the corners, and (iii) renormalize the resultant 
tensors in some proper way. 



1. Simplified one- directional Id method 

Here we assume that the MPO is invariant under trans- 
position of the horizontal indices plus complex conjuga- 
tion, see Fig. (11. a). This symmetry is taken into account 
in order to produce a very efficient algorithm. Thus, the 
TN must be invariant under the index transposition plus 
complex conjugation both in the horizontal (x) and verti- 
cal (y) directions, see Fig.(ll.b). This requirement must 
be satisfied at all steps in the algorithm. In order to 
achieve this, we consider a renormalized TN that is spec- 
ified at every step by one CTM C, a half-column transfer 
matrix Ti, a half-row transfer matrix T2, and the MPO 
tensor see Fig.(ll.b). Notice that tensors Ti and T2 
can also be interpreted as the tensors of two MPSs re- 
spectively in the horizontal and vertical directions (see 
e.g. Fig. (14)). The algorithm implements what we call 
X- and 7/-moves. These are as follows: 
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FIG. 12: (color online) x-move for the simplified one- 
directional Id method, see text. Open indices in reduced 
density matrices are shown in red. 



1. x-move: 

(a) Insertion . We insert one new column in the 
TN, as shown in Fig.(12.a). 

(b) Absorption. We absorb the new column to- 
wards the left, as indicated in Fig.(12.a,b). 
At this step we produce the new (unrenormal- 
ized) tensors C and T2. 

(c) Renormalization . This is done by means of 
the isommetry U as shown in Fig.(12.e). Im- 
portantly, this isommetry is obtained from the 
singular value decomposition ot the CTM C as 
in Fig.(12.c). As a result of this, we obtain the 
renormalized CTM C and half-row transfer 
matrix T2, see Fig.(12.e). Notice that U is, in 
fact, the unitary matrix that diagonalizes the 
reduced density matrix from Fig.(12.d). Fi- 
nally, since the system is symmetric with re- 
spect to horizontal transposition plus complex 
conjugation, we can just use the new tensors 
to obtain the complete TN as shown in the 
diagram of Fig.(ll.b). 

2. y-move: 

(a) Insertion . We insert one new row in the TN, 
as shown in Fig. (13. a). 

(b) Absorption. We absorb the new row towards 
up, as indicated in Fig.(13.a,b). At this tep 
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FIG. 13: (color online) y-uiove for the simplified one- 
directional Id method, see text. Open indices in reduced 
density matrices are shown in red. 



we produce the new (unrenormalized) tensors 
C and fi. 

(c) Renormalization . This is done by means of 
the isommetry V as shown in Fig.(13.e). This 
isommetry is obtained again from the singu- 
lar value decomposition of the CTM C, see 
Fig. (13. c). As in the previous case, V is the 
unitary operator that diagonalizes the reduced 
density matrix in Fig.(13.d). Finally, we use 
the new tensors to obtain once again the com- 
plete TN as in Fig.(ll.b). 

The whole algorithm now proceeds by iterating these 
two moves until convergence of some relevant quantity 
(e.g. the spectrum of singular values of the CTM). This 
approach is valid for computing imaginary time evolu- 
tion (and thus ground states) under an MPO with the 
required symmetry properties, e.g. the quantum Ising or 
antiferromagnetic Heisenberg and XX models described 
in Sec. II. E. The leading number of operations of the 
metod is Oijc')^ where x is the rank of the CTM C. This 
approach has a number of advantages: it is extremely 
efficient in finding ground states (the total number of 
operations is very low), and also preserves the spatial 
symmetries of the evolution operator at every step. Let 
us also stress the fact that the net result of the y-iaove is, 
in fact, equivalent to (i) a rotation of the 2d TN by 7r/2, 
followed by (ii) an x-move, and then followed by (iii) a 
rotation of the TN by — 7r/2. This relation is possible 
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because of the symmetries of the TN in Fig.(ll.b). 

Let us stress another important and intriguing fact 
about this algorithm, which is that there are no explicit 
truncations. That is, the rank x of the CTM is specified 
from the very beggining in the initial TN, and does not 
grow at all throughout the evolution. This is a key differ- 
ence with other methods such as CTMRG and iTEBD, 
where the analogous to this rank grows at every step 
in the algorithm and thus needs to be truncated at every 
step. Here, though, the situation is more subtle. One can 
think of the truncation begin implemented implicitely. 
The initial choice of rank x foi" the CTM defines already 
this implicit truncation, and then x is preserved all along 
the algorithm. Even if this is conceptually strange, it is 
really not a problem for the method: every time we in- 
sert a column or row, the number of indices of the CTM 
proliferates in one direction while it is kept constant in 
the other, which means that the rank of the CTM can 
not grow. It is this property of the algorithm what imple- 
ments an automatic implicit truncation and, thus, we do 
not need to truncate explicitely at any step. However, in 
Sec. III. D we will see that this property no longer holds in 
2d since the indices of the corresponding corner tensors 
always proliferate in, at least, two directions, so explicit 
truncations will always be needed in that case. 



2. Expectation values 

Once the algorithm has converged, one is in position 
of computing expectation values of local operators. For 
this, one needs to consider the contraction of the TN 
together with the operator for the observable that one 
wishes to calculate. The CTM method explained be- 
fore allows to approximate such a calculation. We show 
this for the general case in Fig. (14) (up to normalization 
of the wavefunction^^^) for one-site operators, two-site 
operators between nearest-neighbours, and also for two- 
point correlation functions. As is well known from MPS 
methods, all these calculations can also be done in O(x^) 
operations. 



C. 2d quantum lattice systems 

In this section we discuss how to generalize the ideas 
from the previous section to the 2d case. This general- 
ization is conceptually straightforward, but there are a 
number of significant differences with respect to the Id 
case that are worth stressing. Let us summarize these 
differences: 

(i) Since more spatial dimensions come into play, we 
have more possibilities to effectively expand the rel- 
evant TN. Thus, we discuss four possible corner al- 
gorithms, as opposed to the three discussed in Id. 

(ii) In any of the approaches the indices in the differ- 
ent tensors proliferate at every step in, at least, two 
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FIG. 14: (color online) Different expectation values, up to 
normalization, for a Id system: (a) one-body operator oi, (b) 
two-body operator 02, and (c) two-point correlation function 
with a separation of seven lattice sites. The normalization is 
done by dividing each one of these expressions by the same 
expression where observables are replaced by identity opera- 
tors. 



different directions. Thus, this growth needs to be 
always truncated explicitely, regardless of the ap- 
proach. 

(iii) The spectrums of the corner tensors are no longer 
directly associated to the spectrums of the reduced 
density matrices of the system. Nevertheless, these 
spectrums still carry important information about 
the correlations in the system. 

(iv) The complexity of each approach depends dramat- 
ically on the level of symmetry and chosen renor- 
malization scheme. Unlike in the Id case, where 
all the methods had the same complexity (namely 
0(x^)), in 2d there are differences in several orders 
of magnitude (see Table (I) in Sec. III. E). 

(v) Our algorithms produce separate tensors for the 
'bra' and 'kef parts of the TN, see Fig.(15). Thus, 
the TN is always a positive-definite object by con- 
struction. This was also the case in Id, but was not 
found to be of special relevance there. However, in 
2d this is important since the truncation approach 
differs from the one used in e.g iPEPS and TERG 
algorithms. In those algorithms, the 'bra' and 'kef 
parts of the TN are approximated simultaneously 
by some common tensors, which in practice breaks 
the positivity requirement and may lead to numer- 
ical instabilities (see e.g. the analysis in Ref. ). In 
a way, the approach in this paper is similar to the 
single-layer methods in Ref. . 

Before entering into details of the methods, let us men- 
tion that another alternative for a 2d algorithm based on 
CTMs is actually possible. Namely, one could use some 
Id method with CTMs to compute effective environments 
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FIG. 15: (color online) (a) Symmetries of the PEPO tensor R 
for the simplified one- directional 2d method. In the general 
case (full one-, two- and three-directional 2d methods) only 
the first equality holds, (b) TN structure for the simplified 
one-directional 2d method, (c) General TN structure for the 
full one-, two- and three-directional 2d methods. 



of a 2d iPEPS, and use it in the context of some tensor 
update to simulate time evolution (e.g. 'simplified' or 
'full' updates ' ). This approach is also valid, and has 
already been explored . However, the algorithms that 
we explain in this section are based on a completely dif- 
ferent approach. 

Here, we first present a 'simplified one-directional 2d 
method'. As in the Id case, this simplified method is 
designed for systems such that the tensor defining the 
PEPO is hermitian with respect to the two spatial di- 
rections independently. The implementation of this al- 
gorithm is quite efficient and the steps will be explained 
in detail. In Appendix A.c, we discuss three different 
'full' 2d methods which are valid for PEPO tensors with- 
out the previous symmetry requirement. In increasing 
level of complexity, these approaches are called 'full one- 
directional, full two-directional, and full three-directional 
2d methods'. We only sketch briefly the main idea behind 
them, and from here the interested reader can infer very 
easily their step-by-step implementation. As we shall see 
in the Appendix, these algorithms are less efficient than 
the simplified directional 2d method, and the complexity 
for each one of them is also very different. In fact, some- 
times the complexity may be higher than the calculation 
of expectation values itself, as we will see in Sec. HI. E. 
Still, we believe that it is important (at least from a con- 
ceptual perspective) to be aware of the existence of all 
these possibilities. 

For 2d quantum lattice systems, the CTMs are gener- 
alized to corner tensors (see Fig. (7)), which are tensors 
with at least three indices, one for each direction of the 



lattice. Thus, different singular value decompositions of 
corner tenors may have different ranks. In practice, one 
can work with an effective rank x foi" all the decomposi- 
tions, and this is the relevant parameter for the 2d meth- 
ods. However, it is also possible to work with separate 
ranks, e.g. x the vertical direction and D for the hor- 
izontal directions. This choice of ranks would be analo- 
gous to the iPEPS algorithm, where we have an iPEPS 
of bond dimension D for which an effective environment 
with tensors of bond dimension x^ is computed (in our 
approach, we have that x' ^ ^^^Y roughly, since we 
produce 'bra' and 'kef tensors independently). For sim- 
plicity, in this paper we choose always the same rank 
when explaining the details of the methods (including 
complexity issues). However, we sometimes choose dif- 
ferent ranks for the different directions in the numerical 
calculations of Sec. IV. Let us also remark that, again and 
as in the Id case, the truncations in the vertical indices 
for the full methods are in fact in the entanglement spec- 
trum of a corner of the 2d quantum lattice system, e.g. 
see again Fig. (7). 

As in the Id case, the procedure for these methods is 
always the same: (i) insertion of tensors, (ii) absorption 
of the inserted tensors towards some direction or towards 
some corner, and (iii) renormalization by means of some 
isommetry or rectangular matrix. Also, in this section 
we use the following notation: the spatial (horizontal) 
directions are called x and whereas the temporal (ver- 
tical) direction is called z (see Fig. (15)). Also, x -indices 
refer to indices that connect the tensors in the TN in 
the X direction (and similar definitions apply for y and 
z-indices). 



1. Simplified one- directional 2d method 

Here we assume that the PEPO is invariant under 
transposition of the x and y indices independently plus 
complex conjugation, see Fig. (15. a). As in the Id case, 
this extra symmetry is taken into account in order to 
produce a very efficient algorithm. Thus, the TN must 
be hermitian in the three directions x, y and see 
Fig.(15.b). This requirement will be satisfied at all steps 
in the algorithm. Therefore, the renormalized TN is spec- 
ified at every step by one corner tensor C (analogue to the 
CTM in l(i), three tensors Ti, T2 and T3 (analogue to the 
half-row and half-column transfer matrices in Id), and 
three tensors X, Y and Z (coresponding to some renor- 
malized iPEPS for each one of the three planes yz^ zx 
and xy)^ see Fig.(15.b) The algorithm now implements 
what we call x-, y- and z- moves. As in Id^ these moves 
are equivalent up to rotations of the whole lattice by 7r/2. 
Thus, for simplicity we explain in detail e.g. the x-move, 
and then explain how the y- and z-moves can be related 
to the x-move by rotations. 

1. x-move: 
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FIG. 16: (color online) x-move for the simplfied one- 
directional 2d method, see text. 



(a) Insertion . We insert one new plane of tensors 
in y and z directions of the TN, as shown in 
Fig.(16.a). 

(b) Absorption. We absorb the tensors from the 
new plane towards the left hand side in the 
X direction, as indicated in Fig. (16. a). At 
this step we produce the new (unrenormal- 
ized) tensors C, Ti,T3 and X, see Fig. (16. b). 

(c) Renormalization . This is done by the isom- 
metries Uy^Uz-, Wy and Wz as shown in 
Fig.(16.d). These are found as follows: first, 
we find the isommetries Uy^Uz, Wy and Wz 
from the singular value decompositions of ten- 
sors (7, Ti and T3 that are shown in Fig.(16.b). 
Then, we perform an explicit truncation in 
the X largest singular values respectively of 
all these decompositions, and find the isom- 



y-move = ^ rot + x-move rot 




FIG. 17: (color online) Relation between the x-, y- and z- 
moves in the simplified one-directional 2d method. 



metrics Uy^Uz-, Wy and Wz- With this, we ob- 
tain the renormalized tensors C\T[^T^ and 
X^ see Fig.(16.d). 

2. y-move: rotate the TN by 7r/2 in the xy plane as 
shown in Fig. (17), and do an x-move. Then, rotate 
the TN back to the original position. 

3. z-move: rotate the TN by 7r/2 in the zx plane as 
shown in Fig. (17), and do an x-move. Then, rotate 
the TN back to the original position. 

The algorithm works again by iteration of the above 
three steps until convergence of e.g. the singular value 
spectrums of the corner tensors. As in the Id case, this 
approach is valid for computing ground states by doing 
imaginary time evolution driven by a PEPO with the 
required symmetry properties, e.g. the quantum Ising, 
Heisenberg and XX models described in Sec. II. E. The 
leading number of operations in this method is O(x^) if 
the different contractions are done by making use of the 
geometric structure of the TN (remember the example 
from Fig. (2)). As in Id^ this approach has the advantage 
of being quite efficient, and also keeps the translational 
invariance of the evolution operator at every step. 

Several remarks are in order. First notice that, as 
hinted previously, the indices in the tensors proliferate 
in two different directions at each move, and thus one 
needs explicit truncations in x- But second, and also 
unlike in Id, this time the different truncations are not 
associated to the truncations in any reduced density ma- 
trix of the system. We thus can think of this method as 
an over-simplified algorithm that implements some 'ac- 
ceptable' truncation scheme in a very efficient way. We 
believe that the isommetries found in this way, despite 
not being the best possible, are still good choices as long 
as the amount of entanglement in the system is not too 
large (in a way similar to the simplified update ). In 
the end, whether this approach is useful or not can only 
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be assessed by numerical simulations, and this is what 
we will do in Sec. IV. Notice also that these truncations 
do not correspond to truncations in any entanglement 
spectrum of the system. For this, we should construct 
explicitely the reduced density matrices of the system, 
thus leading to less efficient but probably more accurate 
algorithms. This is precisely what we discuss in the al- 
gorithms of Appendix A.c. 



2. Expectation values 

As in Id, once the algorithm has converged it is pos- 
sible to compute expectation values of local operators. 
This is shown for the general case in Fig. (18) (again up to 
normalization of the wavefunction) for one-site and two- 
site operators between nearest neighbours as well as for 
a two-point function along the x direction. This time the 
required calculations can all be done in 0{x^^) operations 
by choosing the appropriate order in the TN contraction 
of the TN. The fact that the complexity is 0(x^^) means 
that, for the simplfied one-directional 2d method, this is 
actually the bottleneck of the calculation, whereas for the 
rest of the methods explained in Appendix A.c the bot- 
tleneck is the calculation of the reduced density matrices 
in Fig. (30). 



D. Summary of methods and complexities 

In the previous sections we have discussed two differ- 
ent approaches to simulate quantum lattice systems with 
corner transfer matrices and corner tensors. The first ap- 
proach is the one from Ref.*"^: use a (d — 1) -dimensional 
corner method to approximate effective environents in 
the context of an iPEPS algorithm in d dimensions. This 
is the approach that we discussed in Sec. III. D for 3d 
quantum lattice systems. The other approach, which is 
the one that we presented in detail here, is to imple- 
ment directly a corner method in (d + 1) dimensions to 
approximate the contraction involved in the calculation 
of expectation values of local observables for quantum 
lattice systems in d dimensions, assuming an evolution 
driven by some suitable MPO or PEPO. This has been 
done in Sec. III. B and Sec. III. C, as well as Appendix A. 

The complexity in the algorithms presented in 
Sec. III. B, Sec. III. C and Appendix A is summarized in 
Table I. In Id, all the methods that we have studied here 
have the same complexity, which is also the same as in 
iTEBD and iDMRG. However, in 2d the complexity de- 
pends on the chosen method. As expected, 2d methods 
are harder to implement numerically than Id methods. 
Nevertheless, these methods can also be implemented in 
practice within some limitations 




FIG. 18: (color online) Different expectation values, up to 
normalization, for a 2d system: (a) one-body operator oi, (b) 
two-body operator 02, and (c) two-point correlation function 
with a separation of four lattice sites. The normalization is 
done by dividing each one of these expressions by the same 
expression where observables are replaced by identity opera- 
tors. 



Method 


Id 


2d 


Simplified one-directional 


o(x') 


o(x') 


Full one-directional 


o(x') 


o(x") 


Full two-directional 


o{x') 


oix'') 


Full three-directional 




oix'') 


Expectation values 




o(x") 



TABLE I: Complexity of the methods in Sec.III.B, Sec.III.C 
and Appendix A, defined as the leding term in the number 
of operations. The calculation of expectation values is also 
added for comparison. In Id the three complexities are the 
same, but the total running time for each algorithm is dif- 
ferent because of subleading corrections and constant multi- 
plicative terms. 



IV. BENCHMARK: Id AND 2d 

In what follows we present preliminary numerical re- 
sults for some of the algorithms discussed previously. 
Specifically, we have considered the calculation of some of 
the ground state properties of the ferromagnetic quantum 
Ising model in transverse magnetic field h from Eq.(8), 
both in Id and 2d. The properties of this model are well 
known and hence it is useful as a first benchmark 
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FIG. 21: (color online) Two different entanglement spectrums 
FIG. 19: (color online) Entanglement spectrum of half an ^r the 2d ferromagnetic quantum Ising model in a transverse 
infinite chain for the Id ferromagnetic quantum Ising model ^^1^ ^^r different values of h: for a corner (left), and for a 
in transverse field, for different field values h. half-infinite plane (right). 




I 



FIG. 20: (color online) Critical correlation function Cxxil) 
as a funtion of the spin separation / for the Id ferromagnetic 
quantum Ising model in transverse field at criticality. The 
dashed dotted line is the asymptotic behavior in Eq.(16). 



separation distance / as 

a.(0-^+o(^^) (16) 

for some constant a and large In Fig. (20) we plot the 
calculated value of this correlator for ^ = 50 and h = 1. 
We see that our simulation is indeed able to reproduce 
the assymptotic behavior in Eq.(16) for a large number 
of sites even for this small value of x- 



B. 2d 



for our methods. 



A. Id 

In Id the algorithm in Sec. III. B and the ones in Ap- 
pendix A.a,b seem to produce similar results for this 
model. In Fig. (19) we show a calculation of the entangle- 
ment spectrum of half an infinite chain for different val- 
ues of h (focusing on the first 130 spectral values). This 
is computed from the converged spectrum of eigenval- 
ues of the CTM as in Eq.(2). The calculated spectrums 
coincide with remarkable accuracy with the ones in the 
literature (see e.g. Ref.'^^'^^), and were obtained in just a 
few minutes with very modest computational resources. 

We have also computed the two-point correlation func- 
tion Cxx(})i defined as 

where we have substracted the long-distance term m^, 
with rrix the expectation value of (Jx at one site. This 
correlation function can be computed exactly , and at 
criticality {h = 1) it tends to decay algebraically with the 



In 2d we have computed ground state properties of 
the Hamitonian in Eq.(8) by using the simplified one- 
directional 2d method from Sec.III.C.l. In particular, 
we have computed the ground state energy per site eo, 
as well as the magnetizations per site rux and m^, defined 
respectively as the expectation values of the Gx and 
operators at a given site. For completeness and compar- 
ison to the Id case, we have also computed the entangle- 
ment spectrum of a corner of the 2d lattice, as well as of 
half an infinite plane. In our numerical calculations we 
employed different truncation parameters depending on 
the direction. Hence, we used a truncation parameter D 
for the indices of tensor Z in Fig.(15.b) in the xy plane, 
and a parameter x everywhere else. 

In Fig. (21) we plot the entanglement spectrum for dif- 
ferent values of h both for a corner of the infinite plane, 
as well as for half an infinite plane, as computed with 
(D^x) = (4,4). These have been calculated from the 
corner tensors as in e.g. Fig. (7). The finite value of x 
used in our simulations truncates these spectrums in 4 
and 16 values respectively. Our calculations show that 
the entanglement spectrum ffattens as the critical point 
is approached, thus involving a much larger amount of 
entanglement. Interestingly, we also see that long tails 
in the spectrums tend to be produced close to critical- 
ity. When moving slightly away from criticality, we find 
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FIG. 22: (color online) Energy per site (left) and magneti- 
zation vfix per site (right) for the 2d ferromagnetic quantum 
Ising model in transverse field. The calculation is for the sim- 
plified one-directional 2d method with truncation parameters 
(i^,x) = (4,4) 



a crossover region to a regime wehere the spectrums tend 
to decay very quickly. 

In Fig. (22) we show the energy per site eo and the 
magnetization as a function of the magnetic field h 
for (D^x) = (4,4). They follow the usual behavior for 
this 2d quantum system found by other methods (see 
e.g. Fig. (4) in Ref. ) with good accuracy. In Fig. (23) 
we show the behavior of the magnetization m^, which is 
the order parameter, for (D^x) = (4,2), (4,4) and (4,6). 
As expected, close to criticality the order parameter goes 
to zero according to a critical exponent /3 as 



ruz ~ b{hc — hY 



(17) 



for some constant h. The estimated values of P as well as 
our estimations of the critical point for these values of the 
truncation parameters are shown in Fig. (24). Notice that 
for (D, x) = (4, 2) the value of the critical exponent is 
very close to the mean field one /3mf = 1/2, but this gets 
closer to the known Montecarlo result Pmc ^ 0.327 
as higher x is considered. A comparison of the critical 
exponents and critical points as computed by different 
methods is provided in Table II. We see that this ap- 
proach, even if quite simple, already obtains a value for 
the critical exponent that is compatible with the iPEPS 
approach using CTMs from Ref. \ and a value for the 
critical point already better than the Vertical Density 
Matrix Approach (VDMA) from Ref. ' . 

Let us stress that our 2d results are based on the sim- 
plest possible 2d algorithm from this paper. At this 
point, we wish to remind that this algorithm is sort of an 
over-simplified method because the truncating isomme- 
tries are not computed from reduced density matrices at 
all. Thus, much better accuracies in the critical proper- 
ties of the system are expected if instead one implements 
some of the 'full' methods from Appendix A.c. Therefore, 
the numerical results in this paper should be considered 




FIG. 23: (color online) Magnetization rriz per site for the 
2d ferromagnetic quantum Ising model in transverse field. 
The calculation is for the simplified one-directional 2d method 
with truncation parameters {D,x) — (4, 2), (4,4) and (4,6). 
The inset shows the behavior around the critical region. 
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FIG. 24: (color online) Linear fits to the logarithm of rUz 
around the critical region of Fig. (23) to extract the critical 
exponent /3. The different truncation parameters as well as 
the estimated values of /S and the critical point he are indi- 
cated inside the plots. These fits are done at a region suffi- 
ciently close to criticality, but also slightly far from it, in such 
a way that the inherent mean-field effective behavior of the 
TN ansatz for too large correlation lengths is not observed 



only as indicative, and not as representative of the best 
possible performance achievable by corner methods for 2d 
problems. Nevertheless, we feel that it is quite encourag- 
ing that such an over-simplified approach is already able 
to capture the essential properties of the system within 
some accuracy. The numerical performance of some of 
the full 2d methods from Appendix A.c will be consid- 
ered in a future work^^^. 



V. CONCLUSIONS AND FINAL REMARKS 

In this paper we have explored the practical use of 
CTMs and corner tensors to develop tensor network al- 
gorithms for the classical simulation of quantum lattice 
systems of infinite size. At every renormalization step, 
these methods try to minimize the trucation error by 
either (i) keeping the largest singular values of a CTM 
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Method 




he 


Mean Field Theory 


0.5 


4 


Quantum Montecarlo^^*^ 


0.327 


3.044 


i:) = 3VDMA ' 


- 


3.2 


D = 3 MPS + iPEPS'' 


0.332 


3.06 


D = 3 CTM + iPEPS" 


0.328 


3.04 


D = 2 TERG 


0.333 


3.08 


D = 4, X = 4 simp, one-dir. 2d 


0.325 


3.14 



TABLE II: Critical exponent (3 and critical point he of the 
2d quantum Ising model in a transverse magnetic field, as 
computed by different methods. In our case, he has been 
estimated as the point at which rriz becomes roughly 10~^. 



or corner tensor, or (ii) keeping the largest magnitude 
eigenvalues of some reduced density matrix. The second 
scheme maximizes the fidelity between the unrenormal- 
ized and renormalized reduced density matrices. In some 
cases, e.g. Id systems with some extra symmetries, we 
have seen that the truncation in (i) is equivalent to that 
in (ii). 

We have focused mainly on ground states of Id and 
2d systems, although we have also discussed briefly other 
possibilities {Zd systems, periodic boundary conditions, 
and real time evolution). The methods that have been 
proposed preserve the spatial symmetries of the system, 
including invariance under translations. We have bench- 
marked some of this methods by numerically computing 
several ground state properties of the ferromagnetic spin- 
1/2 quantum Ising model in Id and 2d. These numer- 
ics, which should be regarded as preliminary, are already 
quite encouraging. The algorithms of this paper could be 
a possible alternative to other well-established ways to 
compute ground state properties of quantum lattice sys- 
tems in the thermodynamic limit, such as iDMRG and 
iTEBD in Id, and iPEPS and TERG in 2d. The compu- 
tational complexity of the proposed algorithms has also 
been analized, and we have seen that in 2d it depends 
strongly on the simulation scheme, leading to differences 
in several orders of magnitude. 

From a broader perspective, it would be interesting to 
understand the differences between explicit and implicit 
truncations in the rank of the CTMs. In particular, it 
would be good to know whether one type of truncation 
is, by construction, more accurate than the other. In a 
way, this could be similar to the difference between con- 
ventional DMRG and single-site DMRG , where single- 
site DMRG seems to produce more accurate data with 
the same computational resources. 

The methods explored in this paper can also be gen- 
eralized easily to deal with invariance under translations 
every two (or more) lattice sites. In such a case, one 
just needs to choose the tensors accordingly in such a 
way that everything is compatible with the unit cell 
of the lattice. The different renormalizations need also 



to be implemented in accordance with the lattice peri- 
odicity. Finally, let us mention that these algorithms 
can also be used in the context of further tensor net- 
work generalizations in order to study e.g. systems with 
internal symetries and fermionic quantum lattice 
systems^^"^^. We believe that the methods of this pa- 
per will be useful for the practical implementation and 
development of further tensor network algorithms in the 
future. 
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Appendix A: Full Id and 2d methods 

a. Full one- directional Id method 

In this algorithm the only existing symmetry in the 
TN is the vertical transposition of the tensors plus com- 
plex conjugate, and also translation invariance. Thus, 
this time the relevant TN is represented by two CTMs 
Ci and C2, two half-row transfer matrices T2 and T3, and 
one half-column transfer matrix Ti, see Fig.(ll.c). Notice 
that tensors Ti , T2 and T3 can again be interpreted as the 
tensors of three MPS respectively in the horizontal and 
vertical directions (in the vertical direction there is one 
for the left hand side, and one for the right hand side). 
The general idea now is the same as in the previous al- 
gorithm: insert, absorb and renormalize. However, since 
there are less symmetries than in the previous method, 
the renormalization needs to be done in a different way. 
Here are the main steps of the algorithm: 

1. x-move: 

(a) Insertion . We insert two new columns in the 
TN, as shown in Fig. (25. a). 

(b) Absorption. We absorb one new column to- 
wards the left, and one new column towards 
the right, as shown in Fig. (25. a, b). At this 
step we produce four new (unrenormalized) 
tensors Ci , C2 , T2 and T3 . 

(c) Renormalization . This is done by means 
of two isommetries U and V as shown in 
Fig. (25. d). This time, these isommetries are 
computed by calculating the eigenvalue de- 
composition of the reduced density matrices in 
Fig.(25.c). We have that U renormalizes ten- 
sors Ci and T2, whereas V renormalizes C2 
and T3, see Fig. (25. d). It is important that 
U and V are different, since the TN does not 
have extra symmetries in the horizontal direc- 
tion. 

2. y-move: 
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FIG. 26: (color online) y-move for the full one-directional Id 
FIG. 25: (color online) a^-move for the full one-directional Id method, see text. Open indices in reduced density matrices 



method, see text. Open indices in reduced density matrices 
are shown in red 



are shown in red. 



(a) Insertion . We insert one new row, as shown in 
Fig.(26.a). 

(b) Absorption. We absorb the row towards up, as 
in Fig. (26. a, b). At this tep we produce three 
new (unrenormalized) tensors Ci,C2 and Ti. 

(c) Renormalization . This is done as in Fig. (26. d) 
by means of a matrix P. This matrix is com- 
puted by finding the eigenvalue decomposition 
of the reduced density matrix in Fig. (26. c). 
Importantly, this time PP^ is not equal to 
the X X X identity matrix, so we must do 
the renormalization using P and its inverse 
P~^. Notice, however, that the rank of the 
reduced density matrix in Fig. (26. c) is always 
X' Thus, P is a rectangular matrix, and P~^ 
is the pseudoinverse of P. That is, if the sin- 
gular value decomposition of P is given by 
P = FSG^ ^ then its pseudoinverse is defined 



as P-i = G5-ipt. 



The method proceeds again by iterating the above x- 
and y-moves for as long as required until e.g. conver- 
gence of the eigenvalue spectrums of the reduced den- 
sity matrices. Importantly, we can apply this algorithm 
to compute evolutions by MPOs that are not symmet- 
ric. Regarding complexity, this algoithm is again 0{x^)^ 
with X the rank of the reduced density matrices. Also, 
similar conclusions as in the simplified one-directional Id 
method apply for parameter x in this algorithm: it does 
not grow throughout the evolution and is fixed from the 
beginning of the algorithm. Finally, in the case of having 
an MPO invariant under horizontal transposition plus 
complex conjugation, this algorithm is just equivalent to 
the simplfied directional method, although with a less 
efficient implementation. 
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FIG. 27: (color online) Full two-directional Id method: ab- 
sorptions, see text. Open indices in reduced density matrices 
are shown in red. 



b. Full two- directional Id method: CTMRG revisited 

The following algorithm does no longer use a direc- 
tional approach with x- and y-moves. Instead, it uses 
a complete 'radial'-move, expanding both directions of 
the TN at the same time. It is in fact the CTMRG 
method ^^'''^ but adapted to our case (time evolution of 
a Id quantum lattice system driven by an MPO). As in 
the previous section, the system is assumed not to have 
any spatial symmetry in the horizontal direction (apart 
from translation invariance). Thus, the relevant TN is de- 
scribed again by two GTMs Ci and C2, two half-row TM 
T2 and T3, and one half-column TM Ti, as in Fig.(ll.c). 
Again, the half-row and half-column TMs can be inter- 
preted as defining three different MPSs. The details of 
the method are as follows: 

(a) Insertion . We insert two new rows and two new 
columns, as shown in Fig. (27. a). 

(b) Absorption. We absorb the new tensors towards the 
corners by computing the two (unrenormalized) cor- 
ner transfer matrices Ci and C2, half-row transfer 
matrices T2 and T3 and half-column transfer matrix 
Ti as in Fig.(27.b). 




c: c; 

FIG. 28: (color online) Full two-directional Id method: re- 
duced density matrices, see text. Open indices in reduced 
density matrices are shown in red. 

(c) Renormalization . Compute the reduced density ma- 
trices from Fig. (28), and perform their eigenvalue de- 
composition. From these decompositions, obtain the 
two unitary matrices [/, V and the matrix P as well as 
their inverses W and P~^. Truncate explicitely 
in the x largest eigenvalues in magnitude for each 
one of these decompositions, and obtain the matrices 
[/, y, P and their pseudoinverses. Then, renormal- 
ize the CTMs, half-row transfer matrices and half- 
column transfer matrix as shown in Fig. (29). 

These steps are again repeated until convergence of the 
spectrum of the reduced density matrices. Notice that 
this time the rank of the CTM grows at each iteration of 
the method and hence needs to be truncated explicitely in 
X at every step. The leading number of operations for this 
algorithm is again 0{x^). However, the subleading and 
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FIG. 29: (color online) Full two-directional Id method: renor- 
malizations, see text. 



multiplicative corrections are higher than in the previous 
two methods, and hence the total running time is also 
higher. 



c. Full one-, two- and three- directional 2d methods 

From a generic perspective, the 'full' methods in 2d 
have several things in common. First, the PEPO does 
not need to be hermitian in the spatial directions, so we 
will deal with TNs like the one in Fig.(15.c). Second, the 
renormalizations are implemented by tensors found from 
the eigenvalue decomposition of the relevant reduced den- 
sity matrix of the indices that one wishes to truncate. 
Third, truncations in the z direction are implemented 
by isommetries, whereas in the x and y directions one 
uses rectangular matrices and their pseudoinverses. And 
fourth, and unlike in Id, now explicit truncations are al- 
ways needed. 

There are many different ways of puting the above 
ideas in practice, see Fig. (30). For instance, one could 
think of a full one-directional approach, similar to the 
simplified one, but where the tensors that implement 
the truncations are computed from the reduced den- 
sity matrices as in Fig. (30. a). Another approach is 
a full two-directional approach, which keeps the idea 
of the full two-directional method in Id but being re- 
stricted to the planes in the (2 + l)d lattice. In this 
approach, the relevant reduced density matrices are as in 
Fig. (30. b). In fact, one could even think of combining the 
one-directional approach in one direction with the two- 
directional approach in the other two directions, giving 
rise to another possible way of expanding the system. Fi- 
nally, one can think of a full three-directional approach, 
where the reduced density matrices look as in Fig. (30. c). 
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FIG. 30: (color online) Different reduced density matrices 
for two indices that appear at some step in the (a) full 
one-directional, (b) full two-directional and (c) full three- 
directional 2d methods. Reduced density matrices for four 
neighbouring indices from the corners may also be considered 
in the cases (b) and (c), see text. Open indices are shown in 
red. 



This last approach is the most direct generalization of 
the CTMRG algorithm in Ref. ' to higher dimensions. 
However, this is also the most inneficient alternative. No- 
tice also that in Fig. (30) we consider only reduced density 
matrices of two indices. However, it would also be possi- 
ble to consider reduced density matrices of four neigh- 
bouring indices in the same direction within a corner 
in the full two- and three-directional 2d methods (not 
shown). Such an approach would of course be less effi- 
cient, but also more accurate. 

We expect stability and numerical accuracy of these 
methods to increase as the reduced density matrix cap- 
tures more relevant correlations. But the price that 
one has to pay is the growth in complexity which, in 
practice, makes it very difficult to implement numeri- 
cal simulations in some cases. Assuming calculations of 
reduced density matrices for two indices (and thus two- 
index truncations), we see that the full one-directional 
approach has complexity 0(x^^), the full two-directional 
approach 0{^'^\ and the full three-directional approach 
0(x^^). This is to be compared with the simplified one- 



directional method from the previous section, which has 
complexity 0{x^). 

These approaches are expected to produce more accu- 
rate results than the simplified one-directional approach. 
However, in this paper we have just implemented numer- 
ically the simplfied one-directional 2d approach and seen 
that this already produces sensible results, see Sec. IV. 
Thus, we expect these full methods to work even better, 
despite one has to pay a price in the number of operations 
required for the different contractions. The numerical 
exploration of some of these methods will be presented 
elsewhere^*^^. 



Appendix B: Further posibilities 

We now discuss briefly some other possibilities for 
corner-inspired methods, namely, 3d quantum lattice sys- 
tems, periodic boundary conditions, and real time evolu- 
tion. 



d. 3d quantum lattice systems 

3d quantum lattice systems (or equivalently, M clas- 
sical lattice systems) are also of relevance for the study 
of important physical phenomena, e.g. the emergence of 
fermions and gauge bosons in 3d string- net models^^"^, 
or confinement of quarks in Quantum Chromodynamics. 
In principle, it should be possible to generalize all the 
approaches discussed so far in this paper to the 3d case. 
However, from our experience with the 2d case, we expect 
that the complexity of the algorithms will be quite large 
and thus they may become unpractical. Yet, another pos- 
sibility to deal with a 3d quantum lattice system would 
be to follow the same idea as in Ref. for the 2d case, 
namely, to use some 2d method with corner tensors such 
as the one described in Appendix A.c to approximate the 
environment of a 3d iPEPS in the context of a 3d iPEPS 
algorithm. In fact, this would possibly lead to a quite effi- 
cient algorithm for studying 3d quantum lattice systems 
if combined with a simplified update of the tensors^^^. 
For instance, if the full one-directional 2d method is used, 
then we would have a 3d method of complexity 0{x^^)^ 
which is certainly doable. As in Ref.^*^, such an approach 
may break the traslational invariance by one lattice site 
of the lattice, but it would be straightforward to general- 
ize our methods accordingly. The numerical exploration 
of this approach will also be presented elsewhere 

e. Periodic boundary conditions 

Indeed, the case of periodic boundary conditions can 
be considered as well in the context of full algorithms, 
where the renormalizing matrices are computed by eigen- 
value decompositions of the relevant reduced density ma- 
trices. This is so since, in the end, periodic boundary 
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FIG. 31: (color online) Examples reduced density matrices for 
e.g. full one-directional methods with periodic boundary con- 
ditions in (a) Id and (b) 2d. In the 2d case periodic boundary 
conditions are assumed in one spatial direction only, but these 
could also be in both spatial directions. If necessary, further 
periodic boundary conditions in the temporal (vertical) di- 
rection could also be considered. Open indices are shown in 
red. 



conditions just add some extra indices to the TN di- 
agram that needs to be contracted, see e.g. Fig. (31). 
This means that the complexity of the corresponding al- 
gorithm will increase (such an increase in complexity is 
well-known in the DMRG community, see e.g. Ref.^^'^^ 
and also Ref.^^'^^). From the point of view of the con- 
traction of a 2d TN this provides an alternative to the 
method in Ref. . Let us also remark that the case of 
periodic boundary conditions has already been consid- 
ered in the literature in the context of CTMRG, see e.g. 
Ref."^^. Periodic boundary conditions in the temporal di- 
rection may also be useful to compute thermal properties, 
see e.g. Ref.^o^-io^ 



/. Real time evolution: iTEBD and iPEPS reinterpreted 

So far we have discussed the problem of finding ground 
states, which is done by imaginary time evolution. How- 
ever, one may wonder whether the techniques and meth- 
ods explained here are of any use to compute also the real 
time evolution of a quantum lattice system, and hence its 
dynamical properties. We discuss this in what follows. 

In a real time evolution, the MPO or PEPO driving the 
evolution is usually only hermitian with respect to the 
vertical (time) indices. Thus, if one wished to implement 
a time-evolution method based on CTMs and corner ten- 
sors, 'full' methods (or variations thereof) should be con- 
sidered. Moreover, in a real time evolution the total run- 
ning time is finite (unlike in an imaginary time evolution 
for finding ground states). Therefore, we need to deal 
with a 'full' approach such that the time direction can 
grow independently of the spatial directions. The natu- 
ral options for algorithms are thus the full one-directional 



23 





FIG. 32: (color online) (a) The action of an MPO over an 
MPS can be approximated by e.g. three renormalized tensors, 
(b) Reduced density matrix of the bond index of the MPS 
obtained after applying the MPO, in terms of four CTMs. 
Open indices in the reduced density matrix are shown in red. 

approaches from Appendix A.a,b, or in 2d also a combi- 
nation of the full one-directional 2D approach in the time 
direction and the full two-directional 2d approach in the 
space directions (Appendix A.c). 

For concreteness let us focus on the Id case, for which 
some variation of the full one-directional approach from 
Appendix A. a may be used. Based on the above con- 
siderations, one may be tempted to propose the followi- 
ing real-time evolution algorithm: at a given step, (i) 
repeat many x-moves until convergence (effectively mak- 
ing the size of the system infinite), and (ii) perform one 
7/- move. This can be understood as approximating the 
action of an MPO over an MPS by means of three ten- 
sors, see Fig. (32. a). Then, just repeat these operations 
for as many (finite) number of steps as required. 

The above approach is indeed a possibility. However, 
numerical simulations show that its accuracy at long 
times is not as good as the one that can be obtained 
by using e.g. the standard iTEBD method (results not 



shown). This is understandable, since the truncations 
implemented in this approach are more stringent than 
the ones in iTEBD. 

Let us now think of a different approach: we could 
think again on the action of the MPO over the MPS in 
terms of CTMs. Then, the aim is to truncate in the hor- 
izontal indices by means of some tensor. Following the 
spirit of the full algorithms described earlier, this can 
done by computing the reduced density matrix of the in- 
dices by means of the CTMs Ci and C2, see Fig. (32. b). 
Then, we just find its eigenvalue decomposition. With 
this procedure we obtain a matrix P and its inverse P~^. 
A further truncation in the x largest eigenvalues in mag- 
nitude produces the rectangular matrices P and P~^, 
which we can use to renormalize the horizontal indices of 
the tensors. 

Remarkably, this procedure is nothing but a reformula- 
tion of the the standard iTEBD algorithm for non-unitary 
evolutions ' . To see this, notice that iTEBD proceeds 
by (i) finding the canonical form of the resultant MPS 
after the MPO has been applied, and (ii) truncating its 
bond index in the largest x Schmidt coefficients. It is 
easy to prove that for an MPS in canonical form, the re- 
duced density matrix in Fig. (32. b) is already diagonal. In 
our laguage, this diagonalization is actually implemented 
by the matrices P and P~^. Introducing these tensors in 
the TN actually orthonormalizes and truncates the bond 
indices of the evolved MPS. And this is exactly what the 
iTEBD algorithm does. Also, it is not difficult to imagine 
that if we generalize the same procedure to 2d systems, 
then we obtain nothing but a reformulation of the usual 
iPEPS algortihm but using PEPOs for the time evolu- 
tion. In such a case, the relevant reduced density matri- 
ces can not be computed exactly (since they amount to 
the contraction of a 2d TN) and must be approximated 
somehow. Simple approximations lead to efficient algo- 
rithms, whereas more elaborate approximations are less 
efficient but also more accurate. 

Let us remark that the reinterpret ation that we just 
found of the iTEBD algorithm in terms of CTMs is in 
fact similar to the interpretation of DMRG also in terms 
of CTMs, see Ref.^^'^^ 
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The aim of this analogy is simply to make manifest that 
all these algorithms somehow share the same spirit. 
In fact, the largest eigenvalue of the CTM could be degen- 
erate. Yet, this does not spoil the performace of CTMRG 
if only a finite number of eigenvalues is degenerate. Nev- 
ertheless, at criticality all the eigenvalues tend to be de- 
generate, which translates into a critical slow-down of the 
method^^^ 

We adopt here the same notation as in Ref.^°. 
MPOs were first introduced in Ref. 

Here we have used a first-order Suzuki-Trotter decom- 
position, but it is possible to consider higher-order de- 
compositions to reduce the error, e.g. exp(— iif^r) ^ 
exp(-if^^T/2) exp{-H^6r) exp(-i7^^r/2) + O(^r^). 
This analogy is formally valid for the truncations of the 
indices along the spatial directions. 
^^"^ Again, we remark that the actual expectation value is the 
ratio between these expressions, and the same ones but 
where observable operators are replaced by identities. 
See e.g. Ref. . 



